Monte Carlo Study of Doping Change and Disorder Effect 
on Double Exchange Ferromagnetism 
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Phase diagram and critical properties are studied for three-dimensional double exchange model 
with and without quenched disorder. Employing the Monte Carlo method and the systematic analy- 
sis on the finite-size effect, we estimate the Curie temperature and the critical exponent as functions 
of the doping concentration and the strength of the random potential. The Curie temperature well 
scales to the kinetic energy of electrons in the ground state as expected for this kinetics-driven ferro- 
magnetism. The universality class of this transition is described by the short-range Heisenberg fixed 
point. The results are compared with the experimental results in the colossal magnetoresistance 
manganites. 
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I. INTRODUCTION 

Since Zener's pioneering work^ the double exchange 
(DE) model and its extensions have been studied to un- 
derstand the variety of magnetic and transport properties 
including the colossal magnetoresistance (CMR) in per- 
ovskite manganese oxides The original DE model, 
which contains the single-band electrons interacting with 
the localized spins through the Hund's-rule coupling (the 
Hamiltonian will be explicitly given in Sec. Ill A)) , qualita- 
tively explains the stability of the ferromagnetic metallic 
state and the negative magnetoresistance. There, the 
parallel configuration of localized spins leads to the ki- 
netic energy gain of electrons, and vice versa. The effec- 
tive ferromagnetic interaction between localized spins is 
mediated by the motion of electrons, and is called the 
DE interaction. Despite of the long survey for more 
than half century, however, the quantitative aspects of 
the thermodynamics of this model have not been fully 
explored thus far. This is mainly because the system is 
in the strongly-correlated regime in the sense that the 
Hund's-rule coupling is much larger than the bandwidth 
of electrons. Both thermal and spatial fluctuations are 
crucial in the thermodynamics of this itinerant electron 
system, and they are difficult to handle in a controllable 
manner. 

Monte Carlo (MC) calculation is one of the power- 
ful tools to treat such strong correlations properly. In 
the present problem, since the wavefunction of itiner- 
ant electrons is extended, it is crucial to make a sys- 
tematic finite-size scaling analysis for obtaining a reli- 
able result in the thermodynamic limit. Recently, there 
have been some developments in the MC algorithm which 
enable us to handle larger size systems than ever in a 
reasonable computational time&L& The importance of 
spatial fluctuations and finite-size effects has been ex- 



amined by comparing the MC results with the mean- 
field results or the dynamical mean-field approximation 
(DMFA) results^ The reliable estimate of the Curie tem- 
perature Tq has been obtained by the systematic scaling 
analysis p^aLiS The critical exponents have been also esti- 
mated, and it is shown that the exponents are consistent 
with those of the Heisenberg spin model with short-range 
interactions! 8 ^ 10 : 11 ! 12 

These previous MC calculations have been mainly per- 
formed at the doping concentration x — 0.5 (0.5 electrons 
per site on average) , where the kinetic energy of electrons 
is maximum, namely, Tq becomes the highest. In real 
compounds, for instance, in Lai-xSr^MnOa, the ferro- 
magnetic metallic phase is stabilized at 0.15 x ^ 0.6, 
and Tq becomes maximum at x ~ 0.3 and slightly re- 
duces for x > 0.3ii^ This might be due to an instability 
toward the A-type antiferromagnetic state or the charge- 
ordered state near x = 0.5^i5ii£ These instabilities are 
beyond the simple DE model. The x dependence of Tq 
for 0.1 ^ x <; 0.3 in Lai-^Sr^MnOa has been favor- 
ably compared with the DMFA results^ However, the 
recent MC study has revealed the insufficiency of DMFA 
and the importance of spatial fluctuations as mentioned 
above. One of the purposes of this work is to determine 
Tq for wide regions of x by applying the advanced MC 
technique and to clarify the phase diagram of the DE 
model. We will compare the numerical results with the 
experimental results quantitatively. 

Another purpose of this work is to clarify the disor- 
der effect on the critical properties of the ferromagnetic 
transition in the DE model. Disorder suppresses the 
kinetic energy of electrons and reduces Tq. Although 
this reduction has been also studied by DMFA 18 i 19 i 20 i 21 
and the MC calculation for small size clusters^ here we 
give more precise estimates employing the advanced MC 
method and the systematic finite-size scaling. Through 
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the quantitative comparison between Tc and the kinetic 
energy of electrons, we examine how the kinetics governs 
the DE ferromagnetism in the disordered case. We also 
estimate the critical exponent to clarify the universality 
class in the disordered case. These are also motivated by 
the experiments on the ^4-site substitution in AM11O3 at a 
fixed valence x which indicate the relevance of the chem- 
ical disorder of the ^4-site ionsi^iSii2Si2& We will make a 
quantitative comparison between theory and experiment. 

This paper is organized as follows. In SeclHl we intro- 
duce the DE model including the random potential. The 
MC method as well as the details of numerical conditions 
is also described here. In Sec. IIIII we show the numerical 
results. Tc and the critical exponent are estimated by 
the systematic finite-size scaling in both cases with and 
without the random potential. We discuss the numeri- 
cal results in comparison with the experimental results 
in CMR manganites in Sec. IIVI Section is devoted to 
summary. 

II. MODEL AND METHOD 
A. Model 

The DE model considered here consists of itinerant 
electrons in the single band which interact with localized 
spins at each site through the Hund's-rule coupling. The 
Hamiltonian is given in the formi 

Hue = - ^ t ( c L c iff + n - c -) ~ ^ &i • Si, (1) 

where the first term describes the electron hopping be- 
tween the nearest-neighbor sites and the second term rep- 
resents the Hund's-rule coupling between the Pauli spin 
(7 of electrons and the localized spin S (the coupling is 
ferromagnetic, namely, Jh is positive). Additionally, we 
take account of the the random on-site potential as 

H e = ^ j e l c\ a c ia , (2) 

ia 

where we consider the binary distribution, namely, 
takes the value of ±A with equal probability randomly 
in each site. The total Hamiltonian is given by Tt — 
Hde + He- 

In the following, for simplicity, we consider the limit of 
Jh — > oo. In this limit, the electron spin er is completely 
parallel to the localized spin S in each site, and states 
with er antiparallel to S are projected out. This simplifies 
the model to the effective spinless-fermion model in the 
form 2 

n = -J2 *v ^ + h - c -) + E e 4ci, (3) 

where the transfer integral Uj depends on the relative 
angle of localized spins at i and jth site as 

i = t ( cos % cos % + sin % sin . (4) 



Here, the angles 9i and fa are defined as Sf = 
S sin 6i cos fa , Sf = S sin Oi sin fa , and Sf — S cos 0%, 
where S = \Si\ is the magnitude of the localized spin. 
Thus, the transfer integral becomes a complex vari- 
able whose amplitude is proportional to cos(#i — 9j)/2 
and phase is governed by the so-called Berry phase 
exp{-i(fa - fa)}. 

One more simplification we introduce here is to treat 
the localized spins in the classical limit of S — > 00. Then, 
the present model 10 describes the strong interplay be- 
tween the quantum itinerant electrons and the classical 
localized spins. If the configuration of the localized spins 
{Si} is frozen, the problem is simply the free electrons 
in the random magnetic field. However, in the present 
problem, the localized spins are not the fixed external 
magnetic field but internal degrees of freedom of the sys- 
tem. The configurations of {Si} are determined thermo- 
dynamically through the interaction with itinerant elec- 
trons. Namely, thermal equilibrium is achieved to opti- 
mize {Si} which reconciles the kinetics of electrons and 
the entropy of the localized spins. 



B. Method and Numerical conditions 

In the following, we calculate the critical proper- 
ties of the ferromagnetic transition in model © de- 
fined on the three-dimensional cubic lattice. We em- 
ploy the truncated polynomial-expansion MC (TPEMC) 
method which is recently developed by the authors. 7 
This technique is based on the polynomial-expansion MC 
method^ and has advantage in the computational cost 
by introducing effective truncations in the polynomial ex- 
pansion. Readers are referred to Ref. for the details of 
the algorithm. Using this method, we calculate the phys- 
ical quantities in the finite-size clusters of N = L 3 from 
L = 6 to 16. Systematic analyses on the finite-size cor- 
rections are performed for those series. 

We typically perform 10000 MC samplings for mea- 
surements after 1000 steps for thermalization in the ab- 
sence of the disorder. The results are divided into five 
bins to estimate the statistical errors by the variance 
among the bins. In the presence of the disorder, for a 
given configuration of {ei}, we typically perform 1000 
MC samplings for measurements after 1000 steps for ther- 
malization. We repeat this for typically 16 different con- 
figurations to estimate the errors by taking the random 
average among the results for each {£i}. We take the half 
bandwidth W — 6i at Jh = A = as an energy unit. 

We choose twisted boundary conditions to reduce the 
ground state degeneracy in the finite-size systems^ 7 - by 
introducing the magnetic flux <I> = 7r/4, tt/2, and 37r/4 
in the x, y, and z direction, respectively. The magnetic 
flux $ is included by the so-called Peierls factor in the 
phase of the transfer integrals as Uj — > ty exp(i$/L). 
This enables us to reach the converged results in smaller 
order of the polynomial expansion. We confirmed that, 
for the range of the parameters in the present work, the 
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polynomial expansion up to the 8th order is enough to 
obtain the 'exact' results in the limit of the infinite order. 

In the TPEMC calculations, we take the threshold val- 
ues in the truncations as 10 -4 for the matrix product and 
1CP 5 for the trace operation. We confirm that the trun- 
cations with these thresholds do not change the results 
beyond the errors. 

In the following, we calculate the magnetization at a 
fixed doping concentration function of tempera- 

ture, (x = 1 — where the bracket denotes the 

thermal average for the grand canonical ensemble.) For 
that purpose, we need to control the chemical potential 
at each temperature because the band structure changes 
as magnetic correlations develop in this system^ (The 
case of x = 0.5 is special because x is fixed for /i = due 
to the particle-hole symmetry.) Instead of this laborious 
procedure, we fix the chemical potential /i so that the 
target value of x is realized at T ~ Tc. This leads to 
T-dependence of x, however, it is small in the parameter 
range of interests, and does not harm the results because 
the magnetization is not so sensitive to the small change 
of x. We confirmed that the magnetization does not show 
any x dependence beyond its errors for this small devia- 
tion of x. The values of the chemical potential are taken 
as fi = 0.095, 0.205, and 0.310 for x = 0.4, 0.3, and 0.2, 
respectively. The results are shown with the errorbars 
of x which represent the small change of x in the corre- 
sponding temperature range. 



III. RESULTS 
A. Phase diagram in the absence of the disorder 

First, we study the pure case without the random po- 
tential, namely, A = 0. Figure ^ shows the system-size 
extrapolation of the ferromagnetic component of the spin 
structure factor, Sf — J2ij(^i ' divided by the 

system size N. All the data well scale to N~ 2 / 3 , which 
is consistent with the k 2 dependence of the energy cut- 
off of magnons. 2 — From the extrapolated values, we ob- 
tain the magnetization in the thermodynamic limit as 
M = lim^oo y/Sf/N. 

Figure summarizes the temperature dependence of 
the magnetization in the thermodynamic limit. The data 
for x — 0.5 are taken from Ref. llfj. The results are well 
fitted by the scaling form 
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M oc (Tc — Tf 



(5) 



as shown by the gray curves in this figure. From these 
fittings, we estimate the Curie temperature Tc and the 
critical exponent [3 for each value of x. 

The estimates of Tc are summarized in Fig. 0(a). The 
results are shown also for x > 0.5 by using the particle- 
hole symmetry, namely, To (a:) = Tc(l — x). In Fig. 0(b), 
the x dependence of Tc is compared with the kinetic en- 
ergy in the ground state (perfectly polarized state), Ek, 
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FIG. 1: System-size extrapolation of the spin structure factor 
in the pure case for (a) x = 0.4, (b) x — 0.3, and (c) x = 0.2. 
Symbols at N — > oo show the extrapolated values. 



and the functional form of x(l — x) proposed in Ref. l3ll 
The data are normalized by their value at x — 0.5 for 
comparison. As shown in the figure, the MC results al- 
most scale to the kinetic energy of electrons, which con- 
firms that the DE ferromagnetism is governed by the 
kinetics of electrons. From this scaling, we obtain the 
relation T C /\E K /N\ = 0.13 - 0.15. This enables us to es- 
timate Tc approximately from the ground state quantity 
which is easily obtained in this DE system. 

The estimates of /3 by the scaling fit in Fig.0agree with 
the Heisenberg exponent /3 = 0.365 (Ref. 30) within the 
errors for entire range of x. This is consistent with the 
previous MC study at x = 0.5 which has revealed that the 
universality class of this ferromagnetic transition belongs 
to that of the short-range Heisenberg modelAiS 
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FIG. 2: Temperature dependence of the magnetization in the 
thermodynamic limit. Circles, squares, diamonds, and crosses 
represent the data for x = 0.5, 0.4, 0.3. and 0.2, respectively. 
The data for x — 0.5 are from Ref. Iicl The gray curves are 
the scaling fit by eq. JS) . 
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FIG. 3: (a) x dependence of the Curie temperature Tc- (b) 
Comparison of x dependences of Tc (circles), the kinetic en- 
ergy Ek (gray curve), and x(l — x) (dashed curve). All the 
data are normalized by their value at x = 0.5 for comparison. 



B. Disorder effect 

Next, we study the effect of the random on-site poten- 
tial at x = 0.3. Figure0]shows the system-size extrapola- 
tion of the spin structure factor by varying the strength 
of the random potential A. Even in these disordered 
cases, the data well scale to N~ 2 ^ 3 as in the pure case in 
Fig. ^ This is consistent with recent studies for the spin 
excitation spectrum which predict the fe 2 magnon exci- 
tation even in the presence of disorderi22i2£2i The mag- 
netization calculated from the extrapolated values are 
summarized in Fig. [3] The random potential decreases 
the magnetization because it reduces the kinetic energy, 
namely, the DE ferromagnetic interaction. 

We apply the scaling of eq. © in this disordered case 
also, and obtain the estimates of Tc in Fig. (a). Tc 
appears to scale to A 2 in the weak-disorder regime. In 
Fig. (b) , we compare A dependence of Tc with the ki- 
netic energy of electrons in the ground state, and find 
that Tc well scales to the kinetic energy. This indicates 
that, also in the presence of the disorder, the DE ferro- 
magnetism is governed by the kinetics of electrons. 

From the scaling fit in Fig. [5J the critical exponent (3 is 
also estimated. The results are shown in Fig. [7| The ex- 
ponent is consistent with the Heisenberg value (3 = 0.365 
(Ref.HJ even in the presence of the disorder. This indi- 
cates that the disorder is irrelevant and does not change 
the universality class in this transition. 



IV. DISCUSSIONS 



A. Comparison with Lai zSr^MnOs 

Lai-^Sr^MnOs near x — 0.3 has been considered to 
be a canonical DE system in which many aspects of the 
thermodynamics are successfully explained by the DE 
model d alone. 4 9 Here, we compare the estimates of 
Tc obtained in Sec. lIII Al with the experimental values in 
this compound. Comparing the MC estimate T C MC ^ — 

0.0198VF with the experimental value T c ° xp) = 369K 13 
at x — 0.3, we obtain the half bandwidth as W ~ 1.6eV. 
This value is larger than the estimates by the band cal- 
culations (W ~ i e v) however we note that this is 
the bare value in the case of Jh = and that a large Jh 
generally renormalizes the bandwidth. 

Figure |H1 shows the comparison of x dependence of Tc 



between theory and experiment. We normalize T, 



(MC) 



^(oxp) 



well as the kinetic energy of electrons to agree with T c 
at x = 0.3. We note that, in the range of 0.15 ^ x ^ 

0.3, T c cxp ' ) well scales to the MC results and the kinetic 
energy. This agreement has been claimed in the DMFA 
results. 17 For x ^ 0.4, T c cxp ^ is suppressed and shows a 
deviation from this scaling. This might be due to the 
instability toward the A-type antiferromagnetic state or 
the CE-type charge-ordered state observed near x ~ 0.5 
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FIG. 5: Temperature dependence of the magnetization in 
the thermodynamic limit at x = 0.3. Circles, squares, di- 
amonds, crosses, and triangles represent the data for A = 
0.0,0.05,0.1,0.15, and 0.2, respectively. The gray curves are 
the scaling fit by eq. 10. 
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FIG. 4: System-size extrapolation of the spin structure factor 
at x = 0.3 for (a) A = 0.05, (b) A = 0.1, (c) A = 0.15, and 
(d) A = 0.2. Symbols at N — > oo show the extrapolated 
values. 



in many CMR manganitesii2ii£*i& Our data suggest that 
if such instabilities are absent, the Curie temperature can 
become higher up to ~ 420K. 

For the critical exponents, experimental results in this 
compound are still controversial. The estimates are scat- 
tered from the short-range Heisenberg ones to the mean- 
field ones. 38 i 39 i 40 i 41 i 42 i 43 Our results indicate that if the 
DE interaction plays a dominant role in the ferromag- 
netic transition, the Heisenberg universality class should 
be observediSiiSiiiii? Further experimental studies are de- 



1.00 j 




(b) 


rp /rp (A=0) 

A C' 1 C °- 95 






il K /il K 0.90 






0.85 


x=0.3 




n sn 






0.0 0.1 

A 


0.2 



FIG. 6: (a) T c at x = 0.3 plotted as a function of A 2 . The 
line shows the linear ht. (b) Comparison of A dependences 
of Tc (circles) and the kinetic energy Ek (gray curve) . The 
data are normalized by their value at A = for comparison. 



sired. 



Disorder effect in the A-site substitution 



Many experiments have indicated that the disorder 
is important in the ^4-site substitution in ^4Mn03 at a 
fixed doping concentrationiSi 2 ^ 2 ^^ There, the chemi- 
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CMR manganites. 

In the present MC study, we consider only the diag- 
onal disorder of the random on-site potential. In real 
materials, the ^4-site disorder may affect the electronic 
state in other ways, for instance, as the off-diagonal dis- 
order of the random transfer integrals. In the study of 
the spin excitation spectrum, it was shown that vari- 
ous types of disorder bring about the universal results in 
the spectrumi 32 ! 33 ! 34 We speculate that this is the case 
also in the thermodynamics, namely, that Tq scales to 
the kinetic energy and the critical exponents are for the 
Heisenberg universality class irrespective of the types of 
disorder. 

For further substitution which introduces larger differ- 
ence in ionic radius, experimentally, the ferromagnetic 
state is taken over by the charge-ordered state concomi- 
tant with the Jahn- Teller distortioniiii^i& The phase di- 
agram shows the multicritical behavior which indicates a 
strong competition between different phasesi^^SiSSi Dis- 
order effect on the multicritical phenomena has attracted 
much attentions recently^ 1 ^ 2 ^ 3 - The emergence of the 
charge-ordered state suggests that, finally, another ele- 
ment such as the electron-phonon interaction becomes 
important when we approach to the multicritical point. 
However, we consider that, in the regime far from the 
multicritical point, such another element may be less im- 
portant and the disorder plays a primary role. 



FIG. 8: Comparison between the MC results of Tc (cir- 
cles) and the experimental values in Lai-zSr^MnOs (crosses). 
The kinetic energy in the ground state is also shown (gray 
curve). All the data are normalized to agree with each other 
at x — 0.3 for comparison. Triangles represent the antiferro- 
magnetic transition temperature observed in the low x regime. 
Experimental data are from Ref. 



cal disorder from the random distribution of the A ions 
with different ionic radii is considered to disturb the 
electronic state in Mn-O-Mn network electrostatically 
and structurally. For instance, Tq decreases about 30% 
from Lao.7Sro.3Mn03 to Lao.7Cao.3Mn03, and this can- 
not be explained by the pure DE model since the 
estimated change of the transfer integral is only about 
2%M. The MC results in Sec. ITTTBl shows that T c de- 
creases substantially by the disorder. Figure H3 (b) sug- 
gests that the 30% decrease of Tc might be achieved at 
A - 0.4 - 0.5. This corresponds to A ~ 0.6 - 0.8eV if 
we assume W = 1.6eV as in Sec. II V A1 The estimate of 
the disorder strength appears to be consistent with the 
potential fluctuation due to the alloying effect predicted 
by the band calculations^ 

The disorder effect has been also studied in the spin 
excitation spectrum recently^Si^i^ii The results well ex- 
plain the spectral anomalies which show up in the ^4-site 
substituted materialsi^ 4 ^^^ 7 . This and the above ob- 
servation in the decrease of Tc consistently indicate the 
importance of the disorder in the ^4-site substitution in 



V. SUMMARY 

We have studied the phase diagram and the univer- 
sality class of the ferromagnetic transition in the three- 
dimensional double exchange model with and without the 
random potential. The truncated polynomial-expansion 
Monte Carlo method has been employed to calculate 
large size clusters without uncontrolled or biased ap- 
proximations. The Curie temperature and the critical 
exponent have been estimated by applying the system- 
atic finite-size scaling analysis up to 16 3 site clusters. 
For both changes of the doping concentration and the 
strength of the random potential, we found that the 
Curie temperature Tc well scales to the kinetic energy 
of electrons per site in the ground state Ek/N. This 
is a consequence of the fact that the kinetics of elec- 
trons governs the ferromagnetism in this system. From 
this scaling, we have obtained the approximate relation 
T C /\E K /N\ = 0.13 - 0.15, which is useful since the 
ground state quantity is easily calculated in the double 
exchange system. In both cases with and without dis- 
order, estimates of the critical exponent are consistent 
with that of the Heisenberg spin model with short-range 
interaction. This indicates that the ferromagnetic transi- 
tion in the double exchange systems belongs to the short- 
range Heisenberg universality class. We have compared 
the results with the experimental results in Sr doping in 
Lai-^Sr^MnOs and in the ionic radius control by the 
j4-site substitution in j4Mn03. 
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